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Abstract 

We study the mechanical stability of the weakly interacting Bose gas with 
attractive interactions, and construct a unified picture of the collapse valid 
from the low temperature condensed regime to the high temperature classical 
regime. As we show, the non-condensed particles play a crucial role in de- 
termining the region of stability, even providing a mechanism for collapse in 
the non-condensed cloud. Furthermore, we demonstrate that the mechanical 
instability prevents BCS-type "pairing" in the attractive Bose gas. We extend 

our results to describe domain formation in spinor condensates. 
PACS numbers: 03.75.Fi, 05.30.Jp, 64.70.Fx, 64.60.My 
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I. INTRODUCTION 



The trapped Bose gas with attractive interactions is a novel physical system. At high 
densities such clouds are mechanically unstable; however at low densities they can be stabi- 
lized by quantum mechanical and entropic effects. Stability and collapse has been observed 
in clouds of degenerate 7 Li JO and 85 Rb Q. The collapse of the Bose gas shares many 
features of the gravitational collapse of cold interstellar hydrogen, described by the Jeans 
instability ||; related instabilities occur in supercooled vapors. Theoretical studies of the 
attractive Bose gas, typically numerical, have been limited to zero || or very low tempera- 
ture [||,|§ . Here we give a simple analytic description of the region of stability and threshold 
for collapse valid from zero temperature to well above the Bose condensation transition, and 
thus provide a consistent global picture of the instability. As we see, at finite temperature, 
the phase diagram of the Bose gas includes regions where the non-condensed particles play 
a significant role in the collapse 0. 

Our results are summarized in the phase diagram in Fig. [I] which shows three regions: 
normal, Bose condensed, and collapsed. This third region is not readily accessible exper- 
imentally; the system becomes unstable at the boundary (the solid line in the figure). In 
this figure the condensation and collapse lines actually meet at a finite angle. The gross 
features of the phase diagram can be understood qualitatively with dimensional arguments. 
At low temperatures the only stabilizing force is the zero-point motion of the atoms. This 
"quantum pressure" has a characteristic energy per particle of Eq ~ h 2 V~ 2 / 3 /m, where h 
is Planck's constant, V is the volume in which the cloud is confined, and m is the mass of 
an atom. The attractive interactions which drive the collapse are associated with an energy 
U = h 2 a s n/m, where n is the density, and a s (= — 1.45nm for 7 Li and -20 nm for 85 Rb) is 
the s-wave scattering length. Comparison of U with Eq suggests that for low temperatures 
the density n at collapse should be ~ (V 2 ^ 3 a s )^ 1 , independent of temperature. At high tem- 
peratures, the stabilizing force is thermal pressure, P = —TdS/dV, which is characterized 
by the thermal energy E? ~ kgT. Here comparison of U with E? suggests that at high 
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temperatures the density of collapse should be ~ mkBT/h 2 a s , linear in temperature. The 
crossover between the quantum and classical behaviors (n oc T° and nocT) occurs near T c . 
The collapse, as described here, is a phenomena in which the cloud as a whole participates, 
not just the condensate. 

Our basic approach is to identify the collapse with an instability in the lowest energy 
mode of the system (the breathing mode in a spherically symmetric cloud). As the system 
goes from stable to unstable, the frequency of this mode goes from real to complex, passing 
through zero when the instability sets in. The density response function, x(fc,co>), which 
measures the response of the cloud to a probe at wave vector k and frequency u, diverges 
at the resonant frequencies of the cloud. Therefore, by virtue of the vanishing frequency of 
the lowest energy mode, the collapse is characterized by 

X(k = 2n/L,u = 0)- 1 = 0. (1.1) 

Here k = 2n/L is the wave- vector of the unstable mode, whose wavelength L should be of 
order the size of the system. Equation fll-H ) implicitly determines the line of collapse in the 
temperature-density plane. 

To evaluate the response function analytically, we use a local-density approximation, 
replacing the response function x(/c, u) of the inhomogeneous cloud by that of a gas with 
uniform density n. The response of the uniform gas is evaluated at the same frequency 
and wave-vector as for the inhomogeneous system, and n is given by the central density of 
the atomic cloud. The local-density approximation should be valid for temperatures large 
enough that the thermal wavelength, Ay = {2nTi 2 / raksT) 1 / 2 , is much smaller than the size 
of the trap. In all experiments to date, this condition is satisfied, and we treat kAx as 
a small parameter in our calculation. With this approximation, we calculate the line of 
collapse using the well developed theory of x(k,u) of a uniform gas (reviewed in ||). 

In the experiments on 7 Li, the atoms are held in a magnetic trap with a harmonic 
confining potential V(r) = ^muj 2 r 2 , with uo pa 2n x 145s -1 ||. This potential gives the 
cloud a roughly Gaussian density profile (see Sec. [TV)). The temperature T is typically 50 



times the trap energy Too ~ 7nK, so the parameter kAr, ~ yhu/ksT « 0.14, is small. The 
experiments on 85 Rb use softer traps, ^cu ~ 0.6nK, and colder temperatures T rs 15nK, so 
that kA T « 0.2. 

The goal of our approach is to provide a framework for investigating the interplay of 
condensation and collapse. Although our results are not as accurate as can be obtained nu- 
merically (comparing with previous numerical work ||, we find that are results are always 
well within a factor of two of those calculated using more sophisticated models), the con- 
ceptual and computational advantages of working with a uniform geometry far outway any 
loss in accuracy. Due to their simplicity, the arguments used here provide an essential tool 
to choosing which parameter ranges to investigate in future experiments and computations. 

The tools introduced to discuss the collapse of an attractive gas can also be used to 
describe other instabilities in trapped Bose gases. In Section [V] and |V|, we apply these 
methods to the problem of BCS type pairing between bosons, and towards domain formation 
in spinor condensates. 

II. SIMPLE LIMITS 

A. Zero Temperature 

To illustrate our approach we first consider the stability of a zero temperature Bose 
condensate. The excitation spectrum of a uniform gas is ||10|| : 

/ P \ 2 P 

" 2 = + 9 n- = El (2.1) 

\2m J m 

where g = 4ftfi 2 a s /m. In the attractive case, g < 0, all long wavelength modes with 
k 2 /2m < 2\g\n have imaginary frequencies and are unstable. A system of finite size L only 
has modes with k > 2tc/L, and for larger k has an excitation spectrum similar to Eq. fl2.1Q . 
If \g\n < h 2 Tc 2 /mL 2 , the unstable modes are inaccessible and the attractive Bose gas is 
stable. 
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This information is included in the density response function x(k,co) of the dilute zero- 
temperature gas ||, 

X( k > u ) = 7p%i> ( 2 - 2 ) 

where = k 2 /2m. The poles of x are & t the excitation energies, ±Ek- In particular, 
x(k = 27t/L,uj = 0) diverges when \g\n = h 2 TT 2 /mL 2 . 



B. High Temperature (T » T c ) 

We further illustrate our procedure by calculating the stability of an attractive Bose gas 
at temperature much larger than T c , where thermal pressure is the predominant stabilizing 
force. Quantum effects are negligible in this limit, and the line of collapse is simply the 
spinodal line of the classical liquid-gas phase transition JO] , as characterized by Mermin [|12| . 



We neglect finite size effects, and look for an instability in the uniform gas at zero wavevector, 
k = 0, corresponding to finding where x{k = 0, uj = 0) diverges. The susceptibility x(0, 0) = 
—dn/dfi (where /i is the chemical potential), is proportional to the compressability of the 
system, which diverges when the gas becomes unstable. 

At high temperature we work in the Hartree-Fock approximation, where the density is 
given by the self-consistent solution of 

U = I (2 Jo 3 e^~^ - 1' (2 ' 3) 

with Hartree-Fock quasiparticle energies e p = e p + 2gn; here (3 = 1/kgT. In the classical 
limit — 1), n = e^^~ 29n ^ / Ay. The response x(0, 0) has the structure of the random 

phase approximation (RPA), 

X{k,uj) = - — — , 2.4 

1 - 2gXo{k,uJ) 

where Xo(0,0) = — (dn/dfi) e (where the e are held fixed) is the "bare" response. In the 
classical limit Xo(0,0) = —fin. Since Xo(0,0) is negative, the repulsive system (g > 0) is 
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stable. However, for attractive interactions g < 0), the denominator of Eq. Q2.4Q vanishes 
when 2gxo = 1, which in the classical limit occurs when 2\g\n = k^T . 

The above calculation is only valid well above T c . When ~ h 2 /mL 2 , finite size effects 
start to become important, and a more sophisticated approach is needed. If one blindly used 
the above result near T c one would erroneously find that the instability towards collapse 
prevents Bose condensation from occurring. This difficulty can be avoided by working with 
the finite wave- vector response x{k = 2tt/L, uo = 0), to which we now turn. 



III. DENSITY RESPONSE FUNCTION 

The approximation we shall use for x is to treat the response of the gas in the RPA, 
with the simplifying assumption that the bare response of the condensed and non-condensed 
particles is taken to be the response of a non-interacting system. This approach, employed 
by Szepfalusy and Kondor |13[ in studying the critical behavior of collective modes of a Bose 



gas, and later modified by Minguzzi and Tosi ]nj to include exchange, is simple to evaluate 



analytically, and is valid both above and below T c . It generates an excitation spectrum which 
is conserving and gapless ||. At zero temperature it yields the Bogoliubov spectrum, 



Eq. (|2.1|) , and above T c it becomes the standard RPA with exchange. 
The susceptibility in this approximation has the form, 

= 7i *l + nw 2 c m (3-1) 

(1 - tfXoX 1 - 2 9Xo) ~ VXoXo 

where Xo an d Xo are the condensate and non-condensed particle contributions to the response 
of the non-interacting cloud, 

n(h , f d\ f(q-k/2)-f(q + k/2) 
Xa(k,w) = / — -rg — ■ — — , (3.2a) 

J (27T) d LU - {e q+ k/ 2 - Eg-k/2) 

X C o(k,co) = — (3.2b) 

to - e k LU + e k 

Here is the condensate density, the e k = k 2 /2m are the free particle kinetic energies 
as before, and the Bose factors f(k) are given by (e _/3 ( £fc_ ^ — In Appendix A we 
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briefly review the derivation of this response function. At zero temperature Xq = and the 
susceptibility reduces Eq. ( |2.2| ), while above T c , Xo = 0> an d X reduces to Eq. ( |2.4j) . Figure 
|2] shows the class of diagrams summed in this approximation. 

Expanding Xo m the small parameter kA? (see details in Appendix B), we derive for 
T>T C , 

Air 



gX%(k,0) = -2^- 



arctan |e fc /4/i| 1/2 + gy 2 {e?») - |tt//3/j| 1/2 + 0(kA T ) 



kAn 



(3.3) 



where g v {z) = ^/i 1 ' is the polylogarithm function. For chemical potential /x much larger 
in magnitude than /c#T, the system is classical, and Eq. ( |3.3|) reduces to gx\] = —f3gn, as in 



the Hartree-Fock approach, Sec. [IB. Below T c the chemical potential of the non-interacting 



system vanishes and the response functions are: 



9Xo(k,0) = - 1 ^r + 0((kA T ) ), (3.4a) 



g X c (k, 0) = -16*^. (3.4b) 

Using these expressions we calculate the spinodal line separating the stable and unstable 
regions of Fig. [I] by setting k = 2n/L and solving the equation 

X- 1 cc 1 - g(xo + Wo) - 2g 2 XoXo = 0, (3.5) 

which gives the line of collapse as a function of \x and T (for T > T c ) or as a function of 
n and T (for T < T c ). We use the following relations to plot the instability on the n — T 
phase diagram, 

,' A^3/ 2 (en T>T C 
n = (3.6) 

n + A? 3 C(3/2), T<T C . 



Equations (3.4a) and ( |3.4b[ ) indicate that below T c the noncondensate response Xo scales 



as ~ L, while the condensate response Xo scales as k~ 2 ~ L 2 . For realistic parameters, L 
is the largest length in the problem, so that the condensate dominates the instability except 
when no is much smaller than n. Since the condensate is very localized, even a few particles 



in the lowest mode make hq locally much greater than the density of noncondensed particles. 



From the above equations we calculate the maximum stable value of the condensate 
density no- The line of collapse crosses the line of condensation at a temperature T* = 
h 2 /2mkBL\a\. Above this temperature no condensation can occur. For T < T* the collapse 
limits the density of condensed particles to 



(^o)max decreases monotonically with temperature, from the value (^o)max = tt/4Z/ 2 | a| , 
eventually vanishing at T = T* . Using parameters from experiments we find for the 

Rice 7 Li trap, (JV )^ = £ 3 K)mT x ° = 1700 > and T * = 7 - 5 while for the JILA 85Rb tra P> 
(iVo)^x = 120, and T* = 46 nK, The maximum number of particles vs. temperature for 
the two experiments are plotted in Fig. Q these results are consistent with the experiments, 
and agree quite well with numerical mean- field calculations ||. In particular, our curve 
{No) max (T) for lithium has a slope of —1/2.2 nK at T = 0, which lies between the calculated 
slopes of Davis et al. and Houbiers et al. ||. Although (no) max decreases with temperature, 
the non-condensed density n' increases (n' = A~ 3 £(3/2)). Thus the total density at collapse 
need not be monotonic with temperature (cf. the low temperature region of Fig. Qd). 

Future condensate experiments at higher temperatures and densities should be able to 
study the structure in Eq. (|3.7|) , and map out the phase diagram in Fig. [j]. The rubidium 
experiments are performed at temperatures near T*, where the spinodal line intersects T c , 
and in principle should be able to explore the crossover between the quantum mechanical 
and classical behavior of the instability. The lithium experiments are much further away 
from exploring this regime, and in the current geometry, inelastic processes make such an 
investigation impractical [[RJ. Since T* is proportional to 1/L, a softer trap could be used to 
bring this crossover down to lower temperatures where these difficulties are less severe (see 
Fig. ^). More precise numerical studies at higher temperatures are needed to guide these 
experiments. 



In Fig. ^ we show how the line of instability depends on the size of the system, L. 




(3.7) 
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IV. MODELING THE HARMONIC TRAP 



Most experimental and theoretical results are reported in terms of numbers of particles 
instead of density. By appropriately modeling the density distribution of a harmonically 
trapped gas, we can present our conclusions in such a form. Once the interactions are strong 
enough to modify the density distribution significantly, the system undergoes collapse; thus 
we can take the density distribution to be that of non-interacting particles. For ksT ^> too, 
the density profile is well-approximated by 



where V(r) = mtu 2 r 2 /2 is the confining potential, with characteristic length d = (h/mtu) 1 ^ 2 . 
The density of condensed particles at the center of the trap is n . Above T c , n = 0, and 
below T c , /i = 0. Integrating over space, we have 



The instability occurs in the lowest energy mode of the system, the breathing mode, 
whose wave- vector is proportional to 1/d. In a zero temperature non-interacting gas the 
breathing mode has a density profile Sp oc (2r 2 /d 2 — 3) exp (—r 2 /d 2 ), where r is the radial 
coordinate. In momentum space this distribution is peaked at wave- vector k = 2/d. At 
finite temperature thermal pressure increases the radius of the cloud and the wave vector 
of the breathing mode becomes smaller. Since the response of the non-condensed cloud is 
relatively insensitive to the wave- vector, we look for an instability at k = 2/d. 

The resulting phase diagram, Fig. [|, is similar to that in Fig. [1]. The most significant 
difference is that the line of collapse follows the condensation line (on the scale of the 
figure they appear to coincide over a significant temperature range). This behavior can be 
understood by noting that for trapped particles, condensation results in a huge increase in 
the central density of the cloud (a standard diagnostic of BEC) . 




(4.1) 




(4.2) 
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V. PAIRING 



With minor changes the formalism presented here can be used to investigate the insta- 
bility towards forming loosely bound dimers, or "pairs," the Evans-Rashid transition. Such 



an instability occurs in an electron gas at the superconducting BCS transition [18|, and has 
been predicted by Houbiers and Stoof P,|H?| to occur in the trapped alkalis. The pairing is 
signalled by an instability in the T-matrix of the normal phase ||20|| , which plays the role 
that the density response function plays in the collapse. Again, we simulate the finite size of 
the cloud by looking for an instability at k = 2ir/L, as opposed to k = in a bulk sample. 
In analogy to Eq. (|2.4|) , the T-matrix can be written as a ladder sum, 

T(k,u) = - (5.1) 

1- g^(k,uj) 

In this equation, k is the relative momentum of the pair. The instability towards pairing is 
signalled by T — ► oo, when gS — 1. To the same level of approximation as Eq. (|3.2a|) , the 
medium-dependent part of the "pair bubble" S is 

Setting to = 0, and expanding in small k\r, we find 



g E(k,u; = 0) = -4^- 



l7F ■arctan ( |efc/4/i| ' \ + g 1/2 (e^) - |vr//?/i| 1/2 + 0(kA T ) 



kA T I 1+ \e k /4fi\ 1/2 



(5.3) 



Except for the argument of the arctangent, this expression is identical to twice gxo as given 
in Eq. ( |3.3| ). Since arctangent is a monotonic function, and its argument here is smaller 
than in Eq. (|3.3| ), we see that < 2gXo, which implies that the instability towards collapse 
occurs at a lower density than the pairing instability. Thus we conclude that the pairing 
transition does not occur in an attractive Bose gas. Interestingly, in the classical limit, the 
instabilities towards pairing and collapse coincide. 
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VI. DOMAIN FORMATION IN SPINOR CONDENSATES 



The approach used here to discuss the collapse of a gas with attractive interactions also 
describes domain formation in spinor condensates, and gives a qualitative understanding of 



experiments at MIT |21] in which optically trapped 23 Na is placed in a superposition of two 
spin states. Although all interactions in this system are repulsive, the two different spin 
states repel each other more strongly than they repel themselves, resulting in an effective 
attractive interaction. The collapse discussed earlier becomes, in this case, an instability 
towards phase separation and domain formation. The equilibrium domain structure is de- 



scribed in p2| . Here we focus on the formation of metastable domains. 



The ground state of sodium has hyperfme spin F = 1. In the experiments the system is 
prepared so that only the states |1) = \F = l,mp = 1) and |0) = \F = I, nip = 0) enter the 
dynamics. The effective Hamiltonian is then 

H = Jd 3 r V ^ m V ^ + V(r)4l>i + 9 -f 4^^i. (6.1) 

where ipi (i = 0,1) is the particle destruction operator for the state summation over 
repeated indices is assumed. The effective interactions, = 4nh 2 aij/m, are related to the 
scattering amplitudes a F =o and a F =2, corresponding to scattering in the singlet (Fi + F 2 = 0) 
and quintuplet (Fi + F 2 = 2) channel, by 



a = a n = a 01 = a 10 = a F=2 , (6.2a) 
5a = a n - a 00 = (a F=2 - a F=0 )/3. (6.2b) 

Numerically, a = 2.75nm and 5a = 0.19nm. We introduce g = A-Kh 2 a/m and 5g = 
4irh 2 5a/m. In the mean field approximation, the interaction in Eq. ( |6.1| ) becomes a function 
of n m= o and n, the density of particles in the state |0) and the total density, respectively: 

(Hint) = Jd'r [y - , (6.3) 

which shows explicitly the effective attractive interaction. Initially the condensate is static 
with density n = 10 14 cm -3 , and all particles in state 1 1) . A radio-frequency pulse places half 
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the atoms in the |0) state without changing the density profile. Subsequently the two states 
phase separate and form domains from 10 to 50 /mi thick. The trap plays no role here, so 
we can neglect V(r) in Eq. ( |6.1| ) and consider a uniform cloud. 

Linearizing the equations of motion with an equal density of particles in each state, we 
find two branches of excitations corresponding to density and spin waves |33j], 

(i 2 \ 2 j 2 

2^J +3n- + 0(Sg), (6.4a) 

'-s^omn (6.*) 

Since Sg > spin excitations with imaginary frequencies appear. The mode with the largest 
imaginary frequency grows most rapidly, and the width of the domains formed should be 
comparable to the wavelength A of this mode. By minimizing Eq. (|6.4b|) we find A = 



'27r/nSa = 10/rni, in rough agreement with the observed domain size. 
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APPENDIX A: REVIEW OF THE RPA 



Here we give a brief derivation of the response function x(k,uj), Eq. (|3.1|), used in this 
paper. Generically, the response of a gas, \ — Sn/5U, is the direct response to the per- 

12 



turbation xo plus the response to the mean field generated by the disturbed atoms. For a 
normal gas in the Hartree approximation, 5n = Xo$U + Xa9$ n i while including exchange 
gives 5n = Xo$U + xo^9^ n , as in Eq. (]2.4J) . The RPA amounts to making a simple particular 
approximation to the polarization part xo- For our purposes it suffices to take xo = Xo; the 
response of an ideal gas, Eq. (|3.2a|). 



Generalizing the Hartree approximation to the condensed gas simply requires replacing 
Xo with the response of the non-condensed particles Xo P^ us the response of the condensate 
X c , Eq. (pH). 

Including exchange in the degenerate gas requires some work, since exchange only oc- 
curs in interactions involving non-condensed atoms. A simple technique, demonstrated by 



Minguzzi and Tosi JTj] is to look separately at the change in the density of condensed and 
non-condensed atoms, 5uq and 5n. Within Hartree-Fock these changes are related to the 
applied perturbation 5U by 



which gives the relationship 



\XoJ 



SU + 



f xl zxi ^ 1 



V 2x5 2 X s j 



\ 



)n 
5n 



(Al) 



5n 5n 

X = — + — 



Xo + Xo + QXoXo 



5U ■ 5U (l-^)(l-2«S)-VxSXo 
Diagrammatic expressions for these different approximations are shown in Fig. |2]. 



(A2) 



APPENDIX B: ASYMPTOTIC EXPANSIONS 



In this Appendix we derive the kAj- — > asymptotic expansions for the functions Xq an d 
S, defined by Eqs. ( |3.2a|) and (|5.2|) . These expansions are constructed to be valid for all 
/i. We begin by breaking Xo hito two terms, one containing f(q — k/2) and one containing 
f(q + k/2). After shifting q by ±k/2 and integrating out the angular variables, we have 

' {p + k/2 - q)(p - k/2 + qY 



m 



(2TT) 2 k Jo 



dqq /(g) log 



(p+k/2 + q)(p-k/2-q) 



(Bl) 
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with p = muj/k. We extract the important structure by rewriting the logarithm as an 
integral of the form J dx/x. Scaling all lengths by a multiple of the thermal wavelength, we 
find 

$(k,u) = — w dzl(z) (B2) 

7TKA T Jz- 

/°° 3d ! 

dX X *-Bu 1 » ( B3 ) 



where z± = u/2^/e k T ± kA T /A^. The integral has been characterized by Szepfalusy 
and Kondor |13| . In particular, by expanding the distribution function in terms of Matsubara 
frequencies, one arrives at the z — > asymptotic expansion, 



CO oo 

-2tt X>'X^ R [(-^ + 27rw)- y+1)/2 ] , 

j=l u=l 

where dt(z) is the real part of z. Integration of the leading terms gives Eq. ( |3.3| ). 

Following a similar procedure of integrating out the angular variables we write H as 

B ^ = M? T L dzI ^ < B5) 



with z± = (kA/Ay/7r)[iJl + 2uj/ek ± 1]. Using the expansion in Eq. (p4j ) gives Eq. ( |5.3| ). 
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FIG. 1. Phase diagram, density n [cm -3 ] versus temperature T, for a cloud of attractive 
bosons, a) For a cloud of 7 Li with spatial extent L = 3.15^m, scattering length a s = — 1.45nm, 
corresponding to the experiments in Ref. [1], performed at temperatures T near 300 nK. b) For a 
cloud of 85 Rb with scattering length a s = — 20nm and spatial extent L = 3/xm, corresponding to 
the recent experiments in Ref. [2], performed at T near 15 nK. 

Note the logarithmic scales. The solid line separates the unstable (shaded) region from the stable 
region. The dashed line, representing the Bose condensation transition has been continued into 
the collapsed region to illustrate that the two lines intersect. This diagram is drawn for a uniform 
but finite cloud, but can be applied to harmonically trapped gases by taking n to be the central 
density in the trap. 
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a) X = O + O O + O — 0' w + - 

b) O h = O 

c) o HF = o + o + <n> + - 

d) o H = o + — + — 

e) O hf = O + — + — + O + O 

+ <n> + <U> + TD + J0> 

+ ••• 

FIG. 2. Diagrammatic description of the approximation used for the response function. Lines 
with arrows represent propagators of the non-interacting system. Jagged lines represent the contact 
interaction between particles, a) Generically the response function can be expressed in terms of a 
polarization bubble (the bubble with shaded interior), b) In the Hartree approximation (RPA), the 
polarization bubble is taken to be the response of the non-interacting gas. c) In the Hartree-Fock 
approximation, the polarization bubble involves a sum over repeated interactions, d) In a condensed 
system, the polarization bubble for the Hartree approximation must include terms where one of 
the particles is in the condensate, e) The Hartree-Fock approximation for the condensed gas 
requires writing all exchange graphs. This series is most easily summed by the matrix approach in 
Appendix A. 
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FIG. 3. Scaling of the instability threshold with system size. Solid lines show the maximum 
stable density n max for a given temperature. From the top, the system size is L = 0.3, 3, 30 /xm. 
The other parameters are the same as in Fig. |l|a). The dashed line indicates the Bose-Einstein 
condensation transition. Note the three scaling regimes; at low temperature, n max ~ L 2 , near the 
critical temperature n max ~ L 3//2 , and at high temperatures ra max is independent of L (see Eqs. (5) 
and (6)). 
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FIG. 4. Maximum number of condensed particles as a function of temperature, for the a) 7 Li 
and b) 85 Rb experiments. The (&Ay)° terms in Eq. p. 4a have been included in producing these 
plots. 
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FIG. 5. Phase diagram, number N versus temperature T, for harmonically trapped bosons 
with attractive interactions; a) 7 Li, b) 85 Rb. The parameters used correspond to Fig. [Tj. The open 
circle marks where the spinodal line meets with the BEC phase transition. 
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